rm(list=ls())
pacman::p_load(ggthemes, tidyverse, haven, 
               readxl, estimatr, USAboundaries, 
               sf, interflex, sensemakr, psych, 
               sensemakr, devtools, USAboundaries)

rescale01 <- function(x){
  minx <- min(x, na.rm=T)
  maxx <- max(x, na.rm=T)
  return((x-minx)/(maxx - minx))
}

df <- read_csv('datasets/clean_survey1.csv')

# FIGURES AND TABLES

# Cronbach's alpha loading for main DV
psych::alpha(na.omit(cbind(
  df$policy_fuel_efficiency,
  df$policy_smart_grid,
  df$policy_carbon_tax,
  df$policy_ban_sup,
  df$policy_capture_co2,
  df$policy_electric_energy,
  df$policy_fund_meat_altern,
  df$policy_gas_tax)))

# factor analysis Table B1
fa.out <- psych::fa(r = na.omit(cbind(
  df$policy_fuel_efficiency,
  df$policy_smart_grid,
  df$policy_carbon_tax,
  df$policy_ban_sup,
  df$policy_capture_co2,
  df$policy_electric_energy,
  df$policy_fund_meat_altern,
  df$policy_gas_tax)), 
   nfactors = 1, 
   rotate = "varimax") 
summary(fa.out)
fa.out$loadings
sum(fa.out$loadings^2) / 8

# Correlation Table B2
cortab <- round(cor(na.omit(cbind(
  df$policy_fuel_efficiency,
  df$policy_smart_grid,
  df$policy_carbon_tax,
  df$policy_ban_sup,
  df$policy_capture_co2,
  df$policy_electric_energy,
  df$policy_fund_meat_altern,
  df$policy_gas_tax))),2)
rownames(cortab) <- colnames(cortab) <- c('Fuel Efficiency','Smart Grid','Carbon Tax','Plastics','Capture CO2',
                                          'Clean Energy','Meat Alternatives','Gas Tax')
cortab[upper.tri(cortab,diag=T)] <- ""
cortab<-as.data.frame(cortab)
xtable::xtable(cortab)

# Table 2: Main Model ----
lm1 <- lm_robust(policy_scale ~ sea_level_rise,   df, cluster=fips)
lm2 <- lm_robust(policy_scale ~ sea_level_rise +
                   pid7_r + ideo5_c + age + female + college + 
                   income_60_125 + income_over_125 + income_missing +
                   white,   df, cluster=fips)

texreg::texreg(list(lm1,lm2),
               custom.coef.names = c('Intercept','Susceptibility','Party ID (R)', 
                                     'Conservative', 'Age','Female','College',
                                     'Income 60-125k','Income Over 125k','Income Missing',
                                     'White'),
               custom.model.names = c('Support Policy','Support Policy'),
               stars = c(0.1,0.05,0.01), include.ci=F,
               file = 'tables/table2_appendix_table_main_model.tex')


# Are policy attitudes associated with other susceptibilities? 
# Table 4 

summary(m1 <- lm_robust(policy_scale ~ large_fires +
                          pid7_r + ideo5_c + age + female + college + 
                          income_60_125 + income_over_125 + income_missing+ white, 
                        df, cluster=fips))
summary(m2 <- lm_robust(policy_scale ~ heat +
                          pid7_r + ideo5_c + age + female + college + 
                          income_60_125 + income_over_125 +income_missing+ white, 
                        df, cluster=fips))
summary(m3 <- lm_robust(policy_scale ~ wetbulb +
                          pid7_r + ideo5_c + age + female + college + 
                          income_60_125 + income_over_125 +income_missing+ white, 
                        df, cluster=fips))
summary(m4 <- lm_robust(policy_scale ~ farm_crop_yields +
                          pid7_r + ideo5_c + age + female + college + 
                          income_60_125 + income_over_125 +income_missing+ white, 
                        df, cluster=fips))
texreg::texreg(list(m1, m2, m3, m4),
               custom.coef.names = c('Intercept','Susceptibility Fires',
                                     'Party ID (R)', 
                                     'Conservative', 'Age','Female','College',
                                     'Income 60-125k','Income Over 125k',
                                     'Income Missing',
                                     'White','Susceptibility Heat','Susceptibility Wet Bulb',
                                     'Susceptibility Crop Yields'),
               stars = c(0.1,0.05,0.01), include.ci=F,
               custom.model.names = c('Support Policy','Support Policy','Support Policy','Support Policy'),
               reorder.coef=c(1,2,12:14,3:11),
               file = 'tables/table4_otherdisasters.tex')



# contextual controls ----
# Table C3
# various controls for R2
lm1 <- lm_robust(policy_scale ~ sea_level_rise + 
                   pid7_r + ideo5_c + age + female + college + 
                   income_60_125 + income_over_125 + income_missing + 
                   white,   df, cluster=fips)
lm2 <- lm_robust(policy_scale ~ sea_level_rise + pct_white + pct_college + pct_unemployed + median_income + 
                   pid7_r + ideo5_c + age + female + college + 
                   income_60_125 + income_over_125 +income_missing+ 
                   white,   df, cluster=fips)
lm3 <- lm_robust(policy_scale ~ sea_level_rise + pct_white + pct_college + pct_unemployed + median_income + popden +
                   pid7_r + ideo5_c + age + female + college + 
                   income_60_125 + income_over_125 +income_missing+ 
                   white,   df, cluster=fips)
lm4 <- lm_robust(policy_scale ~ sea_level_rise + pct_white + pct_college + pct_unemployed + median_income + popden + total_pop.x +
                   pid7_r + ideo5_c + age + female + college + 
                   income_60_125 + income_over_125 +income_missing+ 
                   white,   df, cluster=fips)
texreg::texreg(list(lm1, lm2, lm3, lm4),
               custom.coef.names = c('Intercept','Susceptibility',
                                     'Party ID (R)', 
                                     'Conservative', 'Age','Female','College',
                                     'Income 60-125k','Income Over 125k','Income Missing',
                                     'White','Pct White','Pct College','Pct Unemp','Median Income',
                                     'Population Density','Total Population'),
               custom.model.names = c('Support Policy','Support Policy','Support Policy','Support Policy'),
               stars = c(0.1,0.05,0.01), include.ci=F,
               file = 'tables/appendix_table_contextualcontrols_main_model.tex')


# robustness checks: sensitivity analysis -----
# Table C4
# sensitivity analysis
model <- lm(policy_scale ~ sea_level_rise +
              pid7_r + ideo5_c + age + female + college + 
              income_60_125 + income_over_125 + income_missing + 
              white,   df)

modelse <- lm_robust(policy_scale ~ sea_level_rise +
                       pid7_r + ideo5_c + age + female + college + 
                       income_60_125 + income_over_125 + income_missing + 
                       white,   df, cluster=fips)

# Appendix Table 
sensitivity <- sensemakr(model=model,
                         treatment = "sea_level_rise",
                         benchmark_covariates = c("pid7_r"),
                         se = sqrt(diag(vcov(modelse))),
                         kd = 1:3)

# short description of results
summary(sensitivity)
ovb_minimal_reporting(sensitivity, format = "latex")

# Robustness Checks: using alternate measure of risk  ----
lmMoody <- lm_robust(policy_scale ~ american_communities_county_risk +
                   pid7_r + ideo5_c + age + female + college + 
                   income_60_125 + income_over_125 + income_missing + 
                     white,  df, cluster=fips)

texreg::texreg(list(lmMoody),
               custom.coef.names = c('Intercept','Susceptibility (Moodys)','Party ID (R)', 
                                     'Conservative', 'Age','Female','College',
                                     'Income 60-125k','Income Over 125k','Income Missing',
                                     'White'),
               custom.model.names = c('Support Policy'),
               stars = c(0.1,0.05,0.01), include.ci=F,
               file = 'tables/appendix_table_main_model_moodys.tex')

# using NOAA alternate measure of SLR
m1 <- lm_robust(policy_scale ~ affected_housing20,   df, cluster=fips)
m2 <- lm_robust(policy_scale ~ affected_housing20 +
            pid7_r + ideo5_c + age + female + college + 
            income_60_125 + income_over_125 + income_missing + 
              white,   df, cluster=fips)

texreg::texreg(list(m1,m2),digits = 3, 
                  include.ci = FALSE,stars = c(0.01,0.05,0.10),
                  custom.coef.names = c('Intercept','Susceptibility (NOAA)','Party ID (R)', 
                                        'Conservative', 'Age','Female','College',
                                        'Income 60-125k','Income Over 125k',
                                        'Income Missing','White'),
                  custom.model.names = c('Support Policy','Support Policy'),
                  file = 'tables/appendix_table_main_model_noaa.tex')

# second, using zip code level controls -----
lm1 <- lm_robust(policy_scale ~ pct_homes_risk,   df, cluster=zip)
lm2 <- lm_robust(policy_scale ~ pct_homes_risk +
                   pid7_r + ideo5_c + age + female + college + 
                   income_60_125 + income_over_125 + income_missing + 
                   white,   df, cluster=zip)


texreg::texreg(list(lm1,lm2),
               custom.coef.names = c('Intercept','Susceptibility (Zip)','Party ID (R)', 
                                     'Conservative', 'Age','Female','College',
                                     'Income 60-125k','Income Over 125k','Income Missing',
                                     'White'),
               custom.model.names = c('Support Policy','Support Policy'),
               stars = c(0.1,0.05,0.01), include.ci=F,
               file = 'tables/appendix_table_main_model_ziprisk.tex')


# For Table C8 run regression in replication_survey2_forsharing
# Replications with Nationscape and CCES in external R file replication_ns_cces_forsharing.R

# Falsification tests ----

lm1 <- lm_robust(bluelm ~ sea_level_rise + 
                   pid7_r + ideo5_c + age + female + college + 
                   income_60_125 + income_over_125 + income_missing + white, 
                 df, cluster=fips)
lm2 <- lm_robust(rr ~ sea_level_rise + 
                   pid7_r + ideo5_c + age + female + college + 
                   income_60_125 + income_over_125 +income_missing+ white, 
                 df, cluster=fips)
lm3 <- lm_robust(affect_republicans ~ sea_level_rise + 
                   pid7_r + ideo5_c + age + female + college + 
                   income_60_125 + income_over_125 +income_missing+ white, 
                 df, cluster=fips)
lm4 <- lm_robust(affect_prius_drivers ~ sea_level_rise + 
                   pid7_r + ideo5_c + age + female + college + 
                   income_60_125 + income_over_125 +income_missing+ white, 
                 df, cluster=fips)
lm5 <- lm_robust(affect_pickup_drivers ~ sea_level_rise + 
                    pid7_r + ideo5_c + age + female + college + 
                    income_60_125 + income_over_125 +income_missing+ white, 
                  df, cluster=fips)

lm6 <- lm_robust(altruism_scale ~ sea_level_rise + 
                    pid7_r + ideo5_c + age + female + college + 
                    income_60_125 + income_over_125 +income_missing+ white, 
                  df, cluster=fips)
lm7 <- lm_robust(future_orientation ~ sea_level_rise + 
                    pid7_r + ideo5_c + age + female + college + 
                    income_60_125 + income_over_125 +income_missing+ white, 
                  df, cluster=fips)
lm8 <- lm_robust(anthropocentrism_scale ~ sea_level_rise + 
                   pid7_r + ideo5_c + age + female + college + 
                   income_60_125 + income_over_125 +income_missing+ white, 
                 df, cluster=fips)
lm9 <- lm_robust(nature_scale ~ sea_level_rise + 
                   pid7_r + ideo5_c + age + female + college + 
                   income_60_125 + income_over_125 +income_missing+ white , 
                 df, cluster=fips)
lm10 <- lm_robust(recreation_scale ~ sea_level_rise + 
                   pid7_r + ideo5_c + age + female + college + 
                   income_60_125 + income_over_125 +income_missing+ white , 
                 df, cluster=fips)

texreg::texreg(list(lm1,lm2,lm3,lm4,lm5),
               custom.coef.names = c('Intercept','Susceptibility','Party ID (R)', 
                                     'Conservative', 'Age','Female','College',
                                     'Income 60-125k','Income Over 125k','Income Missing',
                                     'White'),
               custom.model.names = c('Blue Lives Matter','Racial Resentment','FT Republicans','FT Prius Drivers','FT Pickup Drivers'),
               stars = c(0.1,0.05,0.01), include.ci=F,
               file = 'tables/appendix_falsification1.tex')

texreg::texreg(list(lm6,lm7,lm8,lm9, lm10),
               custom.coef.names = c('Intercept','Susceptibility','Party ID (R)', 
                                     'Conservative', 'Age','Female','College',
                                     'Income 60-125k','Income Over 125k','Income Missing',
                                     'White'),
               custom.model.names = c('Altruism','Future Orientation','Anthropocentric','Nature Scale','Recreation Scale'),
               stars = c(0.1,0.05,0.01), include.ci=F,
               file = 'tables/appendix_falsification2.tex')

# Robustness Test Table C13 Living Near Water
lm1 <- lm_robust(recreation_scale ~ coastal + 
                    pid7_r + ideo5_c + age + female + college + 
                    income_60_125 + income_over_125 + income_missing + white , 
                  df, cluster=fips)
lm2 <- lm_robust(policy_scale ~ pct_water + 
                   pid7_r + ideo5_c + age + female + college + 
                   income_60_125 + income_over_125 +income_missing+ white , 
                 df, cluster=fips)
texreg::texreg(list(lm1, lm2),
                  custom.coef.names = c('Intercept','Coastal','Party ID (R)', 
                                        'Conservative', 'Age','Female','College',
                                        'Income 60-125k','Income Over 125k','Born Again Christian',
                                        'White','Pct Water'),
                  custom.model.names = c('Policy Scale','Policy Scale'),
                  stars = c(0.1,0.05,0.01), include.ci=F,
                  file = 'tables/appendix_coastal_water.tex')


# Robustness Check: Distance from coast Figure C2
dist <- read_csv('datasets/dist_matrix.csv')
dists <- seq(50,1000, by=50)
save <- as.data.frame(matrix(NA, ncol=11, nrow=length(dists)))
for(i in 1:length(dists)){
  shortest_dist <- dist %>%
    filter(mi_to_county < dists[i])
  counties <- c(unique(shortest_dist$county1), unique(shortest_dist$county2))
  ns <- nrow(df[df$fips %in% counties,])
  lm1 <- lm_robust(policy_scale ~ sea_level_rise + 
                     pid7_r + ideo5_c + age + female + college + 
                     income_60_125 + income_over_125 + income_missing + white, 
                   df[df$fips %in% counties,], cluster=fips)
  save[i,] <- lm1 %>% tidy() %>% filter(term == 'sea_level_rise') %>% mutate(dist = dists[i], n=ns) %>% as.data.frame()
  print(i)
}
plot.out <- save %>%
  ggplot(aes(x=V10, y=V2, ymin=V6, ymax=V7)) +
  geom_errorbar(width=0.1) +
  geom_point(shape=21, fill='white',size=3)+ 
  geom_text(aes(x=V10, y=V7 + 0.01,label=V11),size=2) +
  theme_bw() +
  scale_x_continuous(breaks=dists) +
  labs(x='Distance from Coast\n(miles)', y='Coefficient') +
  geom_hline(yintercept=0, color='red',linetype=2) +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=0.75))
plot.out
ggsave(plot.out, width=6, height=4, 
       file='figures/distance_bins.pdf')

# Robustness state fixed effects, region fixed effects, drop one at a time Table C14 ----
lm1 <- lm_robust(policy_scale ~ sea_level_rise  +
                   pid7_r + ideo5_c + age + female + college + income_60_125 + income_over_125 +
                   income_missing + white + factor(region), 
                 df, cluster=fips)
lm2 <- lm_robust(policy_scale ~ sea_level_rise  +
                   pid7_r + ideo5_c + age + female + college + income_60_125 + income_over_125 +
                   income_missing+ white+factor(state_x), 
                 df, cluster=fips)
texreg::texreg(list(lm1, lm2),
               custom.coef.names = c('Intercept','Susceptibility','Party ID (R)', 
                                     'Conservative', 'Age','Female','College',
                                     'Income 60-125k','Income Over 125k','Income Missing',
                                     'White'),
               custom.model.names = c('Policy Scale','Policy Scale'),
               stars = c(0.1,0.05,0.01), include.ci=F, omit.coef = 'factor',
               file = 'tables/appendix_state_region_FEs.tex')


# Drop one state at a time Figure C3
states <- unique(df$state_x)
results_policy <- data.frame(states, estimate = NA,lower=NA, upper=NA)

for(i in 1:length(states)){
  dat <- df %>%
    filter(state_x != states[i])
  results_policy[i,2:4] <- lm_robust(policy_scale ~ sea_level_rise +
                                       pid7_r + ideo5_c + age + female + college + 
                                       income_60_125 + income_over_125 +income_missing+ white,
                                     dat, cluster=fips,alpha = 0.1) %>% tidy()  %>% 
    filter(term == 'sea_level_rise') %>% dplyr::select(estimate,conf.low, conf.high)
  print(i)
}
plot.out <- results_policy %>%
  filter(states != '') %>%
  mutate(states=factor(states, levels=rev(sort(results_policy$states[-3])))) %>%
  ggplot(aes(x=states, y=estimate, ymin=lower, ymax=upper)) +
  geom_errorbar(width=0) +
  geom_point(shape=21, size=1.5, fill='white') + 
  theme_bw() + 
  geom_hline(yintercept=0, linetype=2, color='red') +
  coord_flip() +
  labs(x='',y='Coefficient')
ggsave(plot.out, width=4, height=6, 
       filename = 'figures/appendix_leave_one_out.pdf')


# east coast vs. west coast
west <- lm_robust(policy_scale ~ sea_level_rise +
                    pid7_r + ideo5_c + age + female + college + 
                    income_60_125 + income_over_125 +income_missing+ white, 
                  df %>% filter(state_x %in% c('CA','WA','OR')), 
                  , cluster=fips)
east <- lm_robust(policy_scale ~ sea_level_rise +
                    pid7_r + ideo5_c + age + female + college + 
                    income_60_125 + income_over_125 +income_missing+ white, 
                  df %>% filter(!state_x %in% c('CA','WA','OR')), 
                  , cluster=fips)
df$west_coast <- ifelse(df$state_x %in% c('CA','WA','OR'), 1, 0)
coastal <- lm_robust(policy_scale ~ sea_level_rise*west_coast +
                    pid7_r + ideo5_c + age + female + college + 
                    income_60_125 + income_over_125 +income_missing+ white, 
                  df,  cluster=fips)
texreg::texreg(list(coastal),
               custom.coef.names = c('Intercept','Susceptibility','West Coast','Party ID (R)', 
                                     'Conservative', 'Age','Female','College',
                                     'Income 60-125k','Income Over 125k','Born Again Christian',
                                     'White','SLR * West Coast'),
               custom.model.names = c('Policy Scale'),
               stars = c(0.1,0.05,0.01), include.ci=F,
               file = 'tables/appendix_coastal_interaction.tex')


# Controlling for Hurricanes -----
m1 <- lm_robust(policy_scale ~ sea_level_rise + HRCN_AFREQ +
            pid7_r + ideo5_c + age + female + college + 
            income_60_125 + income_over_125 +income_missing+ white, df, 
          , cluster=fips)

# hurricane 2 -- Cat 3 storms
m2 <- lm_robust(policy_scale ~ sea_level_rise + hurricane_cat3 +
                  pid7_r + ideo5_c + age + female + college + 
                  income_60_125 + income_over_125 +income_missing+ white, df, 
                , cluster=fips)

texreg::texreg(list(m1, m2),
               custom.coef.names = c('Intercept','Susceptibility','Hurricane Incidence','Party ID (R)', 
                                     'Conservative', 'Age','Female','College',
                                     'Income 60-125k','Income Over 125k','Born Again Christian',
                                     'White', 'Hurricane Cat 3'),
               custom.model.names = c('Policy Scale','Policy Scale'),
               reorder.coef = c(1,2,3,13, 4:12),
               stars = c(0.1,0.05,0.01), include.ci=F,
               file = 'tables/appendix_hurricane_control.tex')


### 
# check hurricane disaster declaration -----
# https://www.fema.gov/openfema-data-page/disaster-declarations-summaries-v2
###

m1 <- lm_robust(policy_scale ~ sea_level_rise + n_hurricane_decl_70_24 +
                  pid7_r + ideo5_c + age + female + college + 
                  income_60_125 + income_over_125 +income_missing+ white,   df, cluster=fips)
m2 <- lm_robust(policy_scale ~ sea_level_rise + n_hurricane_decl_10_24 +
                  pid7_r + ideo5_c + age + female + college + 
                  income_60_125 + income_over_125 +income_missing+ white,   df, cluster=fips)

texreg::texreg(list(m1, m2),
               custom.coef.names = c('Intercept','Susceptibility','Hurr Disaster Cnt 70-24','Party ID (R)', 
                                     'Conservative', 'Age','Female','College',
                                     'Income 60-125k','Income Over 125k','Income Missing',
                                     'White', 'Hurr Disaster Cnt 10-24'),
               custom.model.names = c('Policy Scale','Policy Scale'),
               reorder.coef = c(1,2,3,13, 4:12),
               stars = c(0.1,0.05,0.01), include.ci=F,
               file = 'tables/appendix_hurricane_disaster.tex')




